
\documentclass[11pt]{beamer}
\usepackage[UTF8,scheme=plain]{ctex}
\usepackage{listings}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\usetheme{Boadilla}

\usepackage{framed} % 可以用 \begin{shaded}，即背景色块
\definecolor{shadecolor}{rgb}{0.9,0.9,0.9}

\lstset{
	basicstyle=\footnotesize,
    flexiblecolumns,
    showstringspaces = false,
	numbers=left,
    numberstyle=\tiny,
    numbersep=5pt,
	breaklines,%自动换行
	frame=shadowbox,%加边框
	backgroundcolor = \color{red!30!green!30!blue!30},
	%背景色：蓝灰色
	language=C++, % 语言为 c++
	commentstyle=\color{red!80!green!80!blue!80},% 注释用蓝灰色字体标出
	keywordstyle=\color{blue!70}, % 高亮部分用蓝色字体
	tabsize = 4
}

\begin{document}
	\author{ 路毅 \hspace{0.3cm} 曲阜师范大学 }
	\date{\number\year 年 \number\month 月 \number\day 日}
	\title{计算物理第２部分：常用数值方法}
	%\subtitle{}
	%\logo{}
	%\institute{}
	%\date{}
	%\subject{}
	%\setbeamercovered{transparent}
	%\setbeamertemplate{navigation symbols}{}
	\begin{frame}[plain]
	\maketitle
\end{frame}

\kaishu

\begin{frame}
\frametitle{目录}
	\begin{itemize}
		\item {\color{blue} 第四章：数值求导与多项式拟合 }
		\vspace{1cm}
		\item 第五章：数值积分
		\vspace{1cm}
		\item 第六章：方程求根
        \vspace{1cm}
        \item 第七章：寻找函数极值
	\end{itemize}
\end{frame}

\begin{frame}{导数的数值近似}
我们常常需要计算一个函数的导数。
于是我们写出函数的泰勒展开
\begin{eqnarray}
f(a+h) &=& f(a) + hf'(a) + \frac{f''(a)}{2}h^2 + \frac{f'''(a)}{6} h^3 +\cdots \\
f(a-h)  &=& f(a) - hf'(a) + \frac{f''(a)}{2}h^2 - \frac{f'''(a)}{6} h^3 +\cdots
\end{eqnarray}
其中，$h$ 是一个很小的数。
我们可以用如下公式计算 $f'(a)$，
\begin{eqnarray}
	f'(a) = \frac{ f(a+h) - f(a) }{h} + O(h), \\
	f'(a) = \frac{ f(a) - f(a-h) }{h} + O(h).
\end{eqnarray}
所以，我们计算
\begin{eqnarray}
	f'(a) \approx \frac{ f(a+h) - f(a) }{h} , \\
	f'(a) \approx \frac{ f(a) - f(a-h) }{h} .
\end{eqnarray}
作为 $f'(a)$ 的近似，误差为 $O(h)$。
\end{frame}

\begin{frame}{导数的数值近似}
我们也可以设计
\begin{eqnarray}
 f'(a) = \frac{ f(a+h) - f(a-h) }{2h} + O(h^2),
\end{eqnarray}
并用相应的近似公式来计算导数
\begin{eqnarray}
f'(a) \approx \frac{ f(a+h) - f(a-h) }{2h} , 
\end{eqnarray}
这个公式的误差就是 $O(h^2)$，所以这个公式更精确。
\end{frame}

\begin{frame}{导数的数值近似}
如果我们需要两阶导数，我们也可以类似地构造
\begin{eqnarray}
	f''(a) = \frac{ f(a+h) + f(a-h) - 2 f(a) }{h^2} + \frac{ f^{(4)} (a) }{12} h^2 + O(h^4),
	\label{double-derivative}
\end{eqnarray}
用近似公式来计算 $f''(a)$，
\begin{eqnarray}
f''(a) \approx \frac{ f(a+h) + f(a-h) - 2 f(a) }{h^2},
\end{eqnarray}
其误差约为 $ \frac{ f^{(4)} (a) }{12} h^2 $。
\end{frame}

\begin{frame}[fragile]{数值求导：数值误差}
$h$越小，$h^2$越小，导数越精确？\\
$h$越小，$f(a+h)$ 与 $f(a-h)$越接近，即“两个非常相近的数相减”！

下面我们写一个程序，计算 $e^x$ 在 $x=10$ 处的二阶导数，来观察舍入误差的影响。
我们定义相对误差为
\begin{equation}
	\epsilon_r (x)= \frac{ f''_{cal} (x) - f''_{real} (x) } { f''_{real} (x)},
\end{equation}
其中 $f''_{cal} (x), f''_{real} (x)$ 分别为二阶导数的计算值、真实值。
\begin{lstlisting}
#include<iostream>
using namespace std;

#include<cmath>//包括cmath库，后文中就可以调用库中的exp()函数

int main(){

    double h[10];//步长
    for(int i=0;i<10;i++){
        h[i] = pow(0.1,i+1); //　10个步长分别为0.1,0.01,...,1E-10
    }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{数值求导：数值误差}
\lstset{firstnumber=last}
\begin{lstlisting}
    double x = 10;
    double double_derivative;//两阶导数 f''(x) 的计算值
    double relative_error;//计算值的相对误差

    FILE *fp=fopen("double_derivative_ex.txt","w");//定义一个文件流，打开一个 txt 文件，后面利用这个 txt 文件作图
    if(fp==NULL){//如果打开文件失败
        cout<<" failed to open double_derivative_ex.txt"<<endl;
        exit(1);//结束程序，返回异常信号，exit(1)表示异常退出，exit(0)表示正常退出
    }

    for(int i=0;i<10;i++){
        double_derivative = ( exp(x+h[i]) + exp(x-h[i]) - 2*exp(x) )/h[i]/h[i]; //计算两阶导数
        relative_error = double_derivative/exp(x) -1 ; //计算相对误差
        cout<<" h= "<<h[i]<<"  relative error = "<< relative_error <<endl; //输出到屏幕上
        fprintf(fp, "%lf \t %lf \n", log10( h[i] ), log10( fabs(relative_error) ) ); //输出到文件里,用于画图　\t表示一个tab的空白，\n表示换行	
    }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{数值求导：数值误差}
\lstset{firstnumber=last}
\begin{lstlisting}
    fclose(fp);//关闭文件

    return 0;
}
\end{lstlisting}

\begin{itemize}
\item c++中的FILE类型，即文件流，fp是一个文件流指针。
\item fopen　函数打开一个文件，"w"表是写入内容，"r"表示读取内容。
\item 如果 fopen　打开文件失败，会返回NULL值，否则返还指向该文件的文件流指针
\item fscanf或fprintf函数读取或写入内容。
\end{itemize}
\end{frame}

\begin{frame}[fragile]{数值求导：数值误差}
\begin{figure}
	\centering
    \input{double_derivative_ex.tex}
    %\includegraphics[width=0.8\textwidth]{double_derivative_ex.eps}
	\caption{$ \log_{10} | \epsilon_r | \sim \log_{10}(h)$： $e^x$ 在 $x=10$ 处的两阶导数的相对误差与步长的关系}
	\label{double-derivative-relative-error}
\end{figure}
相对误差在 $h \approx 10^{-4}$ 时，取得最小值。怎样理解这个现象？
\end{frame}

\begin{frame}[fragile]{误差分析}
误差分为两部分，一部分是截断误差，一部分是舍入误差。
\begin{equation}
	\epsilon = \epsilon_{truncation} + \epsilon_{round-off},
\end{equation}
截断误差即公式\ref{double-derivative}中被舍去的部分，大约为
\begin{equation}
	\epsilon_{truncation} \approx \frac{ f^{(4)} (a) }{12} h^2,
\end{equation}
而舍入误差，按我们前面“数值误差”一章的分析，
\begin{equation}
	\epsilon_{round-off} \approx \frac{2f (a) }{ h^2} \times 10^{-16},
\end{equation}
\end{frame}

\begin{frame}[fragile]{误差分析}
所以总误差为
\begin{equation}
	\epsilon = \frac{ f^{(4)} (a) }{12} h^2 + \frac{2f (a) }{ h^2} \times 10^{-16},
\end{equation}
将它对 $h$ 求导数，导数为零时误差取得最小值，得到 $ d \epsilon / dh = 0$ 时，
\begin{equation}
	h = 2.2 \times 10^{-4},
\end{equation}
所以 $h$ 取 $2.2 \times 10^{-4}$ 时，计算误差最小。
这与我们在图中看到的现象是一致的！

所以，我们想通过不断减小$h$的值，来得到尽量精确的导数，这个策略是行不通的。

于是人们开发了其他的策略，比如Richardson外推法。
\end{frame}

\begin{frame}{Richardson外推法求导}
我们注意到，用公式
\begin{eqnarray}
 f'(a) = \frac{ f(a+h) - f(a-h) }{2h} + O(h^2),
\end{eqnarray}
求得的数值导数有如下形式，
\begin{equation}
D(h) = \frac{ f(a+h) - f(a-h) }{2h} = D(0) + a_2 h^2 + a_4 h^4 + \cdots,
\end{equation}
其中$D(h)$表示步长选为h时求得的数值导数，而$D(0)$表示导数的理论值。
常数$a_2,a_4,\cdots$与$h$并没有关系。

那么，我们可以定义
\begin{eqnarray}
D^0_0 &=& D(h) = D(0) + a_2 h^2 + a_4 h^4 + \cdots,\\
D^1_0 &=& D(h/2) = D(0) + a_2 (h/2)^2 + a_4 (h/2)^4 + \cdots,\\
&\cdots& \nonumber\\
D^k_0 &=& D(h/2^k) = D(0) + a_2 (h/2^k)^2 + a_4 (h/2^k)^4 + \cdots.
\end{eqnarray}
利用上面的各行，消掉$a_2, a_4, \cdots$就可以得到越来越精确的导数。
\end{frame}

\begin{frame}{总结：数值求导}

\begin{itemize}
\item 数值求导任务，本质上是用差分来完成的。
\item 差分公式的误差由泰勒展开式的余项来确定。
\item 如果直接使用差分公式，要注意舍入误差。
\item 如果精度要求比较高，可以使用 Richardson 外推法。（了解）
\item Richardson 的方法论：消除差分公式中的低阶误差。
\end{itemize}

\end{frame}

\iffalse
\begin{frame}{Richardson外推法求导}
我们可以定义
\begin{eqnarray}
D^k_m = D^{k+1}_{m-1} + \frac{ D^{k+1}_{m-1} - D^k_{m-1} }{ 4^m -1 }, 
\end{eqnarray}
可以利用这个公式，递推出下表中的D值\\
\begin{center}
\begin{tabular}{cccc}
 $D^0_0$ 	&  &  &  \\  
 $D^1_0$	&	$D^0_1$ &  &  \\  
 $D^2_0$	&	$D^1_1$ & $D^0_2$  	&  \\  
 $D^3_0$	&  	$D^2_1$	& $D^1_2$ 	& $D^0_3$  \\  
 $\cdots$	&	$\cdots$  & $\cdots$ 	& $\cdots$ \\ 
\end{tabular} 
\end{center}
上图中，第一列具有$O(h^2)$精度，第二列具有$O(h^4)$精度，第三列具有$O(h^6)$精度。
比如
\begin{eqnarray}
D^0_1 &=& D^1_0 + \frac{ D^1_0 - D^0_0 }{ 4^1 -1 } \nonumber \\
      &=& D(0) - \frac{1}{4} a_4 h^4 - \frac{5}{16} a_6 h^6 + \cdots.
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{Richardson外推法求导}
\begin{lstlisting}
/*
 * Richardson算法：
 * 首先计算
 *      D^k_0 = (f(x+h/2^k) - f(x-h/2^k))/ ( h/2^{k-1} ), k=0,1,...,n
 * 然后根据递推公式
 *      D^k_m = D^{k+1}_{m-1} + \frac{ D^{k+1}_{m-1} - D^k_{m-1} }{ 4^m -1 }
 * 逐步推得　D^0_n 作为返回值，其误差为　O(h^{2(n+1)}
 *
 * 例：
 *
 * D^0_0
 * D^1_0        D^0_1
 * D^2_0        D^1_1   D^0_2
 *
 * 我们定义一个　3x3　的矩阵 deriv，然后逐列计算矩阵元，最后返回　D^0_2, 即 deriv[2][2] 的值。
 * 
 * n 不能是负数，否则程序终止
 */
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Richardson外推法求导}
\lstset{firstnumber=last}
\begin{lstlisting}
double Richardson_deriv(double x, double h, int n, double (*func)(double) ){

        double coef, y;

        double ** deriv = new double * [ n+1 ];//开辟动态内存，得到新的 (n+1)*(n+1) 矩阵
        for(int i=0;i<n+1;i++) deriv[i] = new double [n+1];
        
        for(int i=0;i<=n;i++){ // 赋值第一列 D^k_0, k=0,1,...,n
                coef = pow(0.5,i);// 1/(2^i)
                deriv[i][0] = ( func(x + h*coef ) - func(x - h*coef) ) / (2* h* coef);
        }
        for(int i=1;i<=n;i++){ // 赋值第　1,...,n 列
                for(int j=i;j<=n;j++){ // j 是行数，i是列数
                        deriv[j][i] = deriv[j][i-1] + (deriv[j][i-1] - deriv[j-1][i-1])/ (pow(4.0,i) -1 );//　这一行与递推公式对应
                }
        }
        //cout<<" error is O(h^"<<2*n+2<<")\n";
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Richardson外推法求导}
\lstset{firstnumber=last}
\begin{lstlisting}
        y = deriv[n][n];
        for(int i=0;i<=n;i++)delete [] deriv[i];//释放动态内存
        delete [] deriv;

        return y;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Richardson外推法求导}
当然，上面的代码只定义了一个函数，要运行，还需要加上头文件，以及主函数，并在主函数中调用，比如
\lstset{firstnumber=last}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>//包括cmath库，才能用exp()函数。

double Richardson_deriv(double x, double h, int n, double (*func)(double) );
    
int main(){
    cout<<"y="<<Richardson_deriv(1,0.1,3,exp)<<endl;
    return 0;
}

double Richardson_deriv(...){
    ...函数定义
}
\end{lstlisting}

\end{frame}

\begin{frame}[fragile]{Richardson外推法求导}
用这个办法，求取$exp(x)$在$x=0,1,2,3,4,5$处的一阶导数，得到的结果如表\ref{table:Richardson}，可以看到，仅仅推到$D^0_2$，误差已经小到$10^{-10}$量级。
\begin{table}[h]
\centering
\begin{tabular}{cccccc}
\hline 
$x$ & $D^0_0$ & $D^0_1$ & $D^0_2$ & error of $D^0_2$  \\ 
\hline 
0 & 1.00166750 & 0.99999979 & 1.00000000 & 0.00000000000 \\ 
1 & 2.72281456 & 2.71828126 & 2.71828183 & 0.00000000001   \\ 
2 & 7.40137735 & 7.38905456 & 7.38905610 & 0.00000000002   \\ 
3 & 20.11902956 & 20.08553274 & 20.08553692 & 0.00000000006   \\ 
4 & 54.68919246 & 54.59813866 & 54.59815003 & 0.00000000017   \\ 
5 & 148.66063807  & 148.41312817 & 148.41315910 &  0.00000000046  \\ 
\hline 
\end{tabular} 
\label{table:Richardson}
\caption{数值求导，Richardson算法：exp(x)的一阶导数}
\end{table}
\end{frame}
\fi

\begin{frame}[fragile]{多项式插值}
已知：$(x_0, y_0), (x_1, y_1), (x_2, y_2), \cdots, (x_N, y_N)$一共$N+1$个数据点，其中$x_i$都不相同。\\
假定$P(x)$是一个多项式，令它通过这$N+1$个数据点，然后根据这个多项式推断一个给定点$x_{N+1}$处的函数值
$f(x_{N+1})$：$x_{N+1}$有可能在两个$x_i$之间（即{\bf 插值}）；也有可能在$x_i$覆盖的范围之外（即{\bf 外推}）。

要通过指定的$N+1$个数据点，多项式至少是 $N$ 阶多项式。
\begin{equation}
P(x) = a_0 + a_1 x + \cdots + a_N x^N,
\end{equation}
所以我们需要定下参数$a_0, a_1, \cdots, a_N$。有了这套参数，即可进行插值和外推。

\end{frame}


\begin{frame}[fragile]{插值公式}
拉格朗日给出了插值公式：
\begin{equation}
P_N(x) = \sum^N_{i=0} y_i \prod_{k \neq i} \frac{ x-x_k }{ x_i - x_k },
\end{equation}

可以写成代码，进行练习。
\end{frame}


\begin{frame}[fragile]{插值公式}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>

/*
 * int n:  number of known data points
 * double x: the interpolation/extrapolation point x coordinate
 */
double LagPoly( int n, double * xdata, double * ydata, double x ){

        double y=0, z;
        for(int i=0;i<n;i++){
                z = 1;
                for(int k=0;k<n;k++){
                        if( k != i ) z *= (x - xdata[k] )/ ( xdata[i] - xdata[k] );
                }
                y += ydata[i] * z;
        }
        return y;
}
\end{lstlisting}

\end{frame}

\begin{frame}[fragile]{插值公式}

\begin{lstlisting}
int main(){

        int n;
        cout<<" number of known data points n = "; cin>>n;
        double x;
        cout<<" interpolation/extrapolation position x = "; cin>>x;
        double * xdata = new double [n];
        double * ydata = new double [n];

        ifstream fp("PolyInput.txt");
        for(int i=0;i<n;i++)
                fp >> xdata[i] >> ydata[i];

        cout<<" y(x) = "<< LagPoly( n, xdata, ydata, x )<<endl;

        fp.close();
\end{lstlisting}

\end{frame}

\begin{frame}[fragile]{插值公式}

\begin{lstlisting}
        ofstream out("PolyCurve.txt");
        for(double x = 0.5; x <=4.5; x += 0.01 ){
                out<< x <<"\t" << LagPoly( n, xdata, ydata, x )<<endl;
        }
        out.close();

        return 0;
}
\end{lstlisting}

\end{frame}

\begin{frame}[fragile]{插值公式}

\begin{figure}
\centering
\includegraphics[width=0.7\linewidth]{LagPoly}
\caption{拉格朗日插值公式}
\label{fig:lagpoly}
\end{figure}

\end{frame}


\iffalse
\begin{frame}[fragile]{Nevill外推法}
拉格朗日公式不是一个很高效的算法，也不方便给出误差估计。\\
Nevill给了一个更高效的算法，只需要$N^2$次操作。\\
定义
\begin{equation}
P_0, P_1, P_2, \cdots, P_N,
\end{equation}
为$N+1$个０阶多项式，分别通过第$0,1,\cdots,N$个点。
很显然，０阶多项式就是常数，
\begin{equation}
P_0 = y_0, P_1 = y_1, \cdots, P_N = y_N,
\end{equation}
而
\begin{equation}
P_{01}, P_{12}, P_{23}, \cdots, P_{N-1,N},
\end{equation}
为$N+1$个1阶多项式，分别通过第${0,1},{1,2},{2,3}\cdots,{N-1,N}$个点。
相应地
\begin{equation}
P_{i(i+1)\cdots(i+m)}
\end{equation}
为m阶多项式，分别通过第${i,i+1,i+m}$个点。

显然，我们真正需要的是$P=P_{01\cdots N}$。
\end{frame}

\begin{frame}[fragile]{Nevill外推法}
然而，这些多项式具有递推关系
\begin{equation}
P_{i(i+1)\cdots(i+m)} = \frac{ (x-x_{i+m}) P_{i(i+1)\cdots(i+m-1)} + (x_i -x) P_{(i+1)\cdots(i+m) } }{ x_i - x_{i+m} },
\end{equation}
很容易证明等式右边确实通过第${i,i+1,i+m}$个点。

举$N=3$为例，我们按下图推得$P_{0123}$，

\begin{tabular}{cccc}
$P_0$	&  &  &  \\   
$P_1$	&$P_{01}$&  &  \\   
$P_2$	&$P_{12}$&$P_{012}$ &  \\  
$P_3$	&$P_{23}$&$P_{123}$ & $P_{0123}$ \\ 
\end{tabular} 
\end{frame}

\begin{frame}[fragile]{算法步骤}
\begin{itemize}
    \item	计算$P_0(x),\cdots,P_3(x)$的值。
    \item	用$P_0,P_1$计算$P_{01}(x)$的值，类似地，计算第二列所有多项式的值。
    \item	计算第三列。
    \item 	计算第四列，即得到$P(x)=P_{0123}(x)$的值。
\end{itemize}
\begin{tabular}{cccc}
$P_0$	&  &  &  \\   
$P_1$	&$P_{01}$&  &  \\   
$P_2$	&$P_{12}$&$P_{012}$ &  \\  
$P_3$	&$P_{23}$&$P_{123}$ & $P_{0123}$ \\ 
\end{tabular} 
\end{frame}

\begin{frame}[fragile]{Nevill外推法}
\begin{lstlisting}
/*
 * Nevill 算法
 *
 * 给定 n 个数据点　(x_data[], y_data[])，指定内插/外推点位置　x，Nevill()函数自动进行多项式拟合，并返回 P(x)值
 *
 * 例：
 * P_0
 * P_1  P_{01}
 * P_2  P_{12}  P_{012}
 *
 * 程序定义一个 3x3 的矩阵(虽然上面是一个下三角阵)，然后逐列计算矩阵元的值, 最后返回 P_{012} 的值。
 *
 * 如果数据点的　x　坐标太靠近，误差可能会比较大，所以程序中断。
 *
 */
double Nevill(int n, double *x_data, double *y_data, double x){
        double y;
\end{lstlisting}
\end{frame}        

\begin{frame}[fragile]{Nevill外推法}
\lstset{firstnumber=last}
\begin{lstlisting}
        //检查数据点　x　坐标，如果太靠近就终止程序
        for(int i=0;i<n;i++){
                for(int j=i+1;j<n;j++){
                        if( fabs(x_data[i]-x_data[j]) < 1E-6 ){
                                cout<<" datapoints' x values are too close.\n";
                                exit(1);
                        }
                }
        }

        double ** tempP = new double * [n];//开辟动态内存
        for(int i=0;i<n;i++) tempP[i] = new double [n];

        for(int i=0;i<n;i++) tempP[i][0] = y_data[i];// 第一列：P_0, P_1, \cdots, P_{n-1}
\end{lstlisting}
\end{frame}        

\begin{frame}[fragile]{Nevill外推法}
\lstset{firstnumber=last}
\begin{lstlisting}
        //use the recurrence relations to derive the final polynomial
        // 用递推公式，推出第二列，第三列。。。
        for(int i=1;i<n;i++){
                for(int j=i;j<n;j++){ 
                        tempP[j][i] = (x - x_data[j])*tempP[j-1][i-1] + (x_data[j-i] - x)*tempP[j][i-1];
                        tempP[j][i] /= x_data[j-i] - x_data[j];
                }
        }

        y = tempP[n-1][n-1];

        for(int i=0;i<n;i++) delete [] tempP[i];
        delete [] tempP;//删除动态内存

        return y;
}

\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Nevill外推法}
我们可以写一个主函数，调用上面写的Nevill函数，实现如下功能：
\begin{itemize}
\item 程序自动从一个数据文件中读取数据
\item 程序自动进行多项式拟合
\item 程序自动给出$x_0,\cdots,x_n$区域内的函数值，并输出到文件中，供画图用。
\end{itemize}
\begin{lstlisting}
int main(){

        int n;//数据个数
        double x,y,dy;

        double *x_data;
        double *y_data;

        FILE *fp=fopen("Nevill_input.txt","r");
        if(fp==NULL){
                cout<<"error: failed to open Nevill_input.txt.\n";
                exit(1);
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Nevill外推法}
\lstset{firstnumber=last}
\begin{lstlisting}
        //读一遍输入文件，并将数据行数输出到屏幕上
        n=0;
        while( fscanf(fp, "%*c%*[^\n]") != EOF ){
                n++;
        }
        cout<<"n="<<n<<endl;

        //开辟动态内存，用于存放数据
        x_data = new double [n];
        y_data = new double [n];
        //读取数据
        fp=fopen("Nevill_input.txt","r");
        for(int i=0;i<n;i++){
                fscanf(fp, "%lf %lf", &x_data[i], &y_data[i]);
        }
        fclose(fp);//结束读取
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Nevill外推法}
\lstset{firstnumber=last}
\begin{lstlisting}
        //为了画图，选取x轴范围
        double xmin = x_data[0], xmax = x_data[0];
        for(int i=0;i<n;i++){
                if( xmin > x_data[i] ) xmin = x_data[i];
                if( xmax < x_data[i] ) xmax = x_data[i];
        }
        xmin -= 1;
        xmax += 1;

        int num_grid = 100;//曲线上的点个数
        double step = (xmax - xmin)/num_grid;//曲线上的点的间距

        //计算每个x点的函数值，并输出到文件中
        fp=fopen("Nevill_curve.txt","w");
        if(fp==NULL){
                cout<<"error: failed to open Nevill_curve.txt.\n";
                exit(1);
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{Nevill外推法}
\lstset{firstnumber=last}
\begin{lstlisting}
        for(int i=0;i<num_grid;i++){
                x = xmin + step*i;
                Nevill(n,x_data,y_data,x,&y);
                fprintf(fp,"%lf   %lf\n", x, y);
        }
        fclose(fp);//结束输出

        return 0;
}
\end{lstlisting}
当然，在main函数前面需要包括输入输出库，声明标准库命名空间，以及声明上面定义的Nevill函数，我们省去了这些内容。

如果需要误差估计，可以将$P_{0123...(n-1)}-P_{0123...(n-2)}$作为误差的估计值给出，上面的代码省去了这一步。
\end{frame}

\begin{frame}[fragile]{Nevill外推法}
如果我们把数据文件"Nevill\_input.txt"中的内容设为
\begin{lstlisting}
1   2
2   3
3   5
4   2
\end{lstlisting}
\begin{figure}[h]
\centering
\includegraphics[width=0.6\textwidth]{Nevill}
\label{figure:Nevill}
\caption{多项式拟合与插值：内维尔算法}
\end{figure}
\end{frame}
\fi

\begin{frame}
\frametitle{目录}
	\begin{itemize}
		\item { 第四章：数值求导与多项式拟合 }
		\vspace{1cm}
		\item { \color{blue} 第五章：数值积分 }
		\vspace{1cm}
		\item 第六章：方程求根
        \vspace{1cm}
        \item 第七章：寻找函数极值
	\end{itemize}
\end{frame}

\begin{frame}{数值积分}

{\bf　Newton-Cortes}：在积分区域内选取一系列等距的节点，在节点之间等长的小区间上做近似积分，然后加起来。\\
\vspace{1cm}
{\bf 高斯积分}: 在积分区域内选取一些特定的{\bf 不等距}的节点，然后用这些节点上的函数值，构造整个积分的近似值。\\
\vspace{1cm}
对于震荡比较厉害的积分，应该用第一类，即Newton-Cortes方法。\\
\vspace{1cm}
对于比较平缓的积分，适用高斯积分，高斯积分可以在很少的资源下得到很高的精度。
\vspace{1cm}
\end{frame}

\begin{frame}{Newton-Cortes积分}
如果需要求取积分
\begin{eqnarray}
I_f = \int^{b}_{a} f(x) dx,
\end{eqnarray}
其中 $I$ 是 $Integral$ 的缩写，就是求函数$f(x)$下方的面积，如下图所示。
\begin{figure}[h]
\centering
    \includegraphics[width=0.6\textwidth]{integration_demon.eps}
    \caption{　f(x)积分示意图 }
    \label{figure:integration_demon}
\end{figure}
\end{frame}

\begin{frame}{Newton-Cortes积分}
我们可以将积分区间 $\left[a, b\right]$ 分割为 $n$ 个均等的小区间
\begin{equation}
\left[x_0=a, x_1\right],\cdots,\left[x_{n-1}, x_n=b \right] \nonumber
\end{equation}
把整个积分分为很多小份
\begin{equation}
I_f = \sum^{n-1}_{i=0} \int^{x_{i+1}}_{x_i} f(x) dx.
\end{equation}
\end{frame}

\begin{frame}{中值法}
我们可以将每个区间的积分近似为
\begin{equation}
\int^{x_{i+1}}_{x_i} f(x) dx \approx h f(\frac{x_i + x_{i+1}}{2}),
\end{equation}
当区间足够小的时候，小区间中的函数值都很接近，所以可以做上面的近似。
总积分则近似为
\begin{eqnarray}
I_f \approx \sum^{n-1}_{i=0} h f(\frac{x_i + x_{i+1}}{2}).
\end{eqnarray}
这个方法取每个小区间的中点的函数值，进行近似，所以叫做中值法。
\end{frame}


\begin{frame}[fragile]{中值法}
\begin{lstlisting}
/*
 * integ_midpoint: 中值法积分
 * 给定区间[a,b]，分割为n份，被积函数为func，返回积分值
 */
double integ_midpoint(double a, double b, int n, double (*func)(double)){

        int i;
        double y=0;
        double h=(b-a)/n;
        for(i=0;i<n;i++){
                y += func( a + (i+0.5)*h );
        }
        y *= h;
        return y;
}

double f(double x){
        return x*x;
}
\end{lstlisting}
\end{frame}

\begin{frame}{梯形法}
我们也可以换一个方法做近似
\begin{equation}
\int^{x_{i+1}}_{x_i} f(x) dx \approx h (f(x_i) + f(x_{i+1}))/2,
\end{equation}
上式右边正是梯形面积公式。
这个近似的几何意义是：用一个梯形的面积对小区间内的积分做近似，所以这个方法叫做梯形法。
总的积分值为
\begin{eqnarray}
I_f &\approx& \sum^{n-1}_{i=0} h (f(x_i) + f(x_{i+1}))/2 \nonumber\\
    &=& h \{ \frac{f(x_0)}{2} + f(x_1) + \cdots + f(x_{n-1}) + \frac{f(x_n)}{2} \}.
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{梯形法}
\begin{lstlisting}
/*
 * integ_trapzoid: 梯形法积分
 * 给定区间[a,b]，分割为n份，被积函数为func，返回积分值
 */
double integ_trapzoid(double a, double b, int n, double (*func)(double)){

        int i;
        double y=func(a)/2;
        double h=(b-a)/n;
        for(i=1;i<=n-1;i++){
                y += func( a + i*h );
        }
        y += func(b)/2;
        y *= h;
        return y;
}
\end{lstlisting}
\end{frame}

\begin{frame}{Simpson方法}
在小区间内，中值法是用一个常数作为函数的近似
\begin{equation}
f(x) \approx a,
\end{equation}
而梯形法是用１阶多项式作为函数的近似
\begin{equation}
f(x) \approx a + bx.
\end{equation}
那么，如果我们想要更精确的近似，可以取三个点
\begin{eqnarray}
&&(x_i, f(x_i)), \nonumber\\
&&(x_i+h/2, f(x_i+h/2)), \nonumber\\
&&(x_i+h, f(x_i+h))，
\end{eqnarray}
进行２阶多项式拟合。
\end{frame}

\begin{frame}{辛普生方法}
利用前文中多项式拟合的拉格朗日公式，可以得到
\begin{eqnarray}
f(x) &\approx& f(x_i) \frac{ (x-x_{h/2})(x-x_h) }{ (-h/2)(-h) } \nonumber\\
    && + f(x_i + h/2) \frac{ (x-x_i)(x-x_h) }{ (h/2) (-h/2)} \nonumber\\
    && + f(x_i + h) \frac{ (x-x_i)(x-x_{h/2}) }{ h \cdot h/2},
\end{eqnarray}
为了方便，我们用$x_{h/2},x_h$代表$x_i + h/2, x_i + h$。
在区间 $\left[x_i, x_{i+1}\right]$ 内积分以后得到
\begin{equation}
\int^{x_i +h}_{x_i} f(x) dx \approx \frac{h}{6}[ f(x_i) + 4f(x_i + h/2) + f(x_i + h) ].
\end{equation}
\end{frame}

\begin{frame}{辛普生方法}
所以有
\begin{eqnarray}
I_f &\approx& \sum^{n-1}_{i=0} \frac{h}{6} \left[ f(x_i) + 4 f(\frac{x_i + x_{i+1}}{2}) + f(x_{i+1}) \right] \nonumber\\
&=& \frac{h}{6} \left\{ f(x_0) + 4f(x_0+\frac{h}{2}) + 2f(x_1) + 4f(x_1+\frac{h}{2}) + \cdots \right. \nonumber\\
 && \left. + 4f(x_{n-1}+\frac{h}{2}) + f(x_n)    \right\}
\end{eqnarray}
这就是Simpson方法。
\end{frame}

\begin{frame}[fragile]
\begin{lstlisting}
/*
 * Simpson：辛普生积分
 * 给定积分区间[a,b]，区间分割为n份，被积函数为func，返回积分值
 */
double Simpson(double a, double b, int n, double (*func)(double)){

        double y1,y2,y,h;
        int i;
        h = (b-a)/n;//步长
        y1 = 0;
        for(i=0;i<n;i++){
                y1 += func(a+ (i+0.5)*h);
        }
        y2 = 0;
        for(i=1;i<n;i++){
                y2 += func(a+ i*h);
        }
        y = 4*y1 + 2*y2 + func(a) + func(b);
        y *= h/6;
        return y;
}
\end{lstlisting}
\end{frame}

\begin{frame}{误差分析}
利用泰勒展开，很容易推导出Newton-Cortes积分公式的误差。
为了方便叙述，我们记$m = x_i + h/2$。
我们写出$f(x)$在$m$附近的
\begin{eqnarray}
f(x) =&& f(m) + f'(m)(x-m) + \frac{f''(m)}{2} (x-m)^2  \nonumber\\ 
        &&+ \frac{ f^{(3)}(m)}{6}(x-m)^3 + \frac{ f^{(4)}(m)}{4!}(x-m)^4 + \cdots
\end{eqnarray}
其中 $f^{(n)}(m)$ 表示函数在 $x=m$ 处的 $n$ 阶导数。
所以有
\begin{eqnarray}
f(x_{i}) = && f(m) - f'(m)\frac{h}{2} + \frac{f''(m)}{2} \frac{h^2}{4}  \nonumber\\ 
        && - \frac{ f^{(3)}(m)}{6}\frac{h^3}{8} + \frac{ f^{(4)}(m)}{4!}\frac{h^4}{16} + \cdots \\
f(x_{i+1}) = && f(m) + f'(m)\frac{h}{2} + \frac{f''(m)}{2} \frac{h^2}{4}  \nonumber\\ 
        && + \frac{ f^{(3)}(m)}{6}\frac{h^3}{8} + \frac{ f^{(4)}(m)}{4!}\frac{h^4}{16} + \cdots 
\end{eqnarray}
\end{frame}

\begin{frame}{误差分析}
$f(x)$在区间$\left[x_i, x_{i+1}\right]$上的积分为
\begin{eqnarray}
\int^{x_{i+1}}_{x_i} f(x)dx &=& f(m)h  + \frac{f''(m)}{24} h^3 + \frac{ f^{(4)}(m)}{1920}h^5 + \cdots 
\label{equation:IntegralOn1Step}
\end{eqnarray}

将$f(x_i), f(m), f(x_{i+1})$的泰勒展开分别带入Newton-Cortes公式，得到它们的展开形式。
中值法
\begin{eqnarray}
\int^{x_{i+1}}_{x_i} f(x)dx &\approx& f(m)h,
\end{eqnarray}
梯形法
\begin{eqnarray}
\int^{x_{i+1}}_{x_i} f(x)dx &\approx& f(m)h + f''(m)h^3/8 + f^{(4)}(m)h^5/(4!2^4) + \cdots , \nonumber\\
\end{eqnarray}
\end{frame}

\begin{frame}{误差分析}
辛普生方法
\begin{eqnarray}
\int^{x_{i+1}}_{x_i} f(x)dx &\approx& f(m)h + f''(m)h^3/24 + f^{(4)}(m)h^5/(1152) + \cdots , \nonumber\\
\end{eqnarray}
中值法、梯形法的误差均为$O(h^3)$，辛普生方法的误差为$O(h^5)$。
\vspace{1cm}

但这是每一步长上的积分误差，所以总误差为：中值法、梯形法$O(h^2)$，辛普生方法$O(h^4)$。
\end{frame}

\begin{frame}{自适应算法}
数$f(x)$震荡的特别厉害的时候，我们需要比较小的步长$h$; \\
如果被积函数$f(x)$非常平缓，较大的步长也可以得到不错的精度。

于是我们希望，程序根据函数的性质，自动选择合适的步长，进行积分。

以辛普生方法为例，上面推导出，它在$[x_i,x_{i+1}]$的误差
\begin{equation}
S^{[x_i,x_i+h]}_f - I^{[x_i,x_i+h]}_f = c h^5,
\end{equation}
其中$S_f$表示新普生方法得到的数值积分，c为常数(我们假定$f^{(4)}(x)$在$[x_i,x_{i+1}]$区间内变化很小。)

那么，如果我们将步长减半，则会得到
\begin{eqnarray}
&& S^{[x_i,x_i+h/2]}_f - I^{[x_i,x_i+h/2]}_f = c h^5/32, \nonumber\\
&& S^{[x_i+h/2,x_i+h]}_f - I^{[x_i+h/2,x_i+h]}_f = c h^5/32,
\end{eqnarray}
即
\begin{eqnarray}
S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - I^{[x_i,x_i+h]}_f = c h^5/16,
\end{eqnarray}
\end{frame}

\begin{frame}{自适应算法}
所以
\begin{eqnarray}
S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - S^{[x_i,x_i+h]}_f = - \frac{15}{16} c h^5,
\end{eqnarray}
即
\begin{eqnarray}
&& S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - S^{[x_i,x_i+h]}_f \nonumber\\
=&& - 15 (S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - I^{[x_i,x_i+h]}_f),
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{自适应算法}
\begin{eqnarray}
&& S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - S^{[x_i,x_i+h]}_f \nonumber\\
=&& - 15 (S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - I^{[x_i,x_i+h]}_f),
\end{eqnarray}
所以，如果我们要求$[x_i,x_{i+1}]$内的积分误差不超过$\epsilon$，
\begin{eqnarray}
| S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - I^{[x_i,x_i+h]}_f | \leq \epsilon
\end{eqnarray}
则只需要求
\begin{eqnarray}
|S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - S^{[x_i,x_i+h]}_f | < 15 \epsilon
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{自适应算法}
我们可以写一个递归函数，实现自适应积分，代码如下
\begin{lstlisting}
/*
 * function Simpson applies the Simson's rule to calculate integration of func() over [a,b], breaking [a,b] into n small pieces.
 */
double Simpson(double a, double b, int n, double (*func)(double)){
        double h = (b-a)/n;

        double y = func(a) + func(b);
        for(int i=0;i<n;i++){
                y += 4*func(a+ (i+0.5)*h);
        }
        for(int i=1;i<n;i++){
                y += 2*func(a+ i*h);
        }
        y *= h/6;
        return y;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{自适应算法}
\lstset{firstnumber=last}
\begin{lstlisting}
int count=0;
/*
 * SelfAdapt_Simpson uses self adapted algorithm to calculate integration of func() over [a,b].
 * It automatically breaks [a,b] into unequal pieces, according to the local situaion of func().
 * the error is controled to be less than epsilon
 */
double SelfAdapt_Simpson(double a, double b, double (*func)(double), double epsilon){

        count++;

        double y1 = Simpson(a, (a+b)/2, 1, func);
        double y2 = Simpson((a+b)/2, b, 1, func);
        double y3 = Simpson(a, b, 1, func);
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{自适应算法}
\lstset{firstnumber=last}
\begin{lstlisting}
        if( fabs(y1+y2-y3) < 15*epsilon ){
                return y1+y2;
        }
        else{
                return SelfAdapt_Simpson(a, (a+b)/2, func, epsilon/2)
                        +
                        SelfAdapt_Simpson((a+b)/2, b, func, epsilon/2);
        }
}
\end{lstlisting}
其中，Simpson函数实现Simpson积分，而SelfAdapt\_Simpson函数实现自适应的Simpson积分。
\end{frame}

\begin{frame}[fragile]{自适应算法}
我们可以做一个例子，比如做如下积分
\begin{equation}
\int^{0.1}_{10} \sin(1/x^2) dx,
\end{equation}
如图所示，被积函数$\sin(1/x^2)$在图左侧振荡得非常厉害，在右侧却非常平缓。
所以在左侧需要很小的步长，在右侧稍大的步长即可。
\begin{figure}[h]
\centering
\includegraphics[width=0.6\textwidth]{integration_demon_self_adapt.eps}
%\caption{自适应辛普生积分：$\int^{10}_{0.1} \sin(1/x^2) dx$}
\label{figure:integration_demon_self_adapt}
\end{figure}
\end{frame}

\begin{frame}[fragile]{自适应算法}
我们分别用两种方式做这个积分，一种是变步长自适应辛普生，另一种是同步长辛普生方法，得到的结果如下表。
在同样的精度要求下，自适应积分只需要十分之一的资源。
\begin{table}[h]
\centering
\begin{tabular}{cccc}
\hline\hline
\multicolumn{2}{c}{自适应} & \multicolumn{2}{c}{普通} \\
\hline 
步数 & 误差 & 步数 & 误差 \\ 
\hline 
1E4 & 4.71E-10 & 1E4 & -2.5E-6 \\ 
\hline 
 &  & 1E5 & -2.24E-10  \\
\hline 
\end{tabular} 
\caption{自适应辛普生积分：$\int^{10}_{0.1} \sin(1/x^2) dx$}
\label{table:self-adapt}
\end{table}
\end{frame}

\begin{frame}{高斯积分：多项式作为近似}
如果在$[a,b]$内用N阶多项式$P(x)$作为$f(x)$的近似，让$P(x)$经过$(x0=a,y0=f(a)), \cdots, (x_N=b, y_N=f(b))$一共N+1个点，则在$[a,b]$区间内任意一处，有
\begin{equation}
f(x) = P(x) + \frac{f^{(N+1)}(\epsilon)}{(N+1)!} (x-x_0)\cdots(x-x_N),
\end{equation}
积分以后误差为
\begin{eqnarray}
\epsilon_P &=& \int^b_a P(x)dx - \int^b_a f(x)dx        \nonumber\\
&=& \int^b_a \frac{f^{(N+1)}(\epsilon)}{(N+1)!} (x-x_0)\cdots(x-x_N) dx \nonumber\\
&\leq& \frac{f^{(N+1)}(\epsilon)}{(N+1)!} (b-a)^{N+2}   
\end{eqnarray}
\end{frame}

\begin{frame}{高斯积分：多项式作为近似}
\begin{eqnarray}
\epsilon_P &=& \int^b_a P(x)dx - \int^b_a f(x)dx        \nonumber\\
&=& \int^b_a \frac{f^{(N+1)}(\epsilon)}{(N+1)!} (x-x_0)\cdots(x-x_N) dx \nonumber\\
&\leq& \frac{f^{(N+1)}(\epsilon)}{(N+1)!} (b-a)^{N+2}   
\end{eqnarray}
$\frac{f^{(N+1)}(\epsilon)}{(N+1)!} (b-a)^{N+2}$随着$N$的增大，会收敛为0。

所以，一般来说，多项式的阶数越高，误差越小。
\end{frame}

\begin{frame}{勒让德多项式}
勒让德多项式是一类微分方程的解
\begin{equation}
C(1-x^2)P - m^2_l P + (1-x^2)\frac{d}{dx}((1-x^2)\frac{dP}{dx}) = 0.
\end{equation}
其中$m_l=0$时，解是勒让德多项式；$m_l\neq0$时，解是连带勒让德多项式。

勒让德多项式为
\begin{equation}
L_k (x) = \frac{1}{2^k k!} \frac{d^k}{dx^k} (x^2-1)^k, k=0,1,2,\cdots
\end{equation}
它具有正交性
\begin{equation}
\int^1_{-1} L_i(x) L_j(x) dx = \frac{2 \delta_{ij} }{2i+1},
\end{equation}
\end{frame}

\begin{frame}{勒让德多项式}
前面几项勒让德多项式为
\begin{eqnarray}
& L_0(x) = 1, \nonumber\\
& L_1(x) = x, \nonumber\\
& L_2(x) = \frac{1}{2}(3x^2-1), \nonumber\\
& \cdots \nonumber
\end{eqnarray}
\end{frame}

\begin{frame}{分解高阶多项式}
前面我们解释了，如果想要提高积分精度，可以用高阶的多项式作为函数的近似。
如果积分区域为$[-1,1]$，要求
\begin{equation}
\int^1_{-1} f(x) dx,
\end{equation}
我们想用$P_{2N-1}$，即2N-1阶多项式来作为$f(x)$的近似，我们可以将它分解为
\begin{eqnarray}
P_{2N-1}(x) = L_{N}(x) P_{N-1}(x) + Q_{N-1}(x),
\end{eqnarray}
这是可以做到的。
左侧有$2N$个参数，右侧$P_{N-1}(x),Q_{N-1}(x)$分别有$N$个参数，这个分解是唯一确定的。
\end{frame}

\begin{frame}{分解高阶多项式}
我们注意到，$L_0(x), \cdots, L_{N-1}(x)$分别为$0,1,\cdots,N-1$阶多项式，所以$P_{N-1}(x)$一定可以表示为
\begin{equation}
P_{N-1}(x) = \sum^{N-1}_{i=0} a_i L_i(x),
\end{equation}
所以
\begin{equation}
\int^1_{-1} L_{N}(x) P_{N-1}(x) dx
= \sum^{N-1}_{i=0} a_i L_N (x) L_i(x) dx = 0,
\end{equation}
\end{frame}

\begin{frame}{高斯积分}
所以有
\begin{eqnarray}
\int^1_{-1} f(x)dx \approx \int^1_{-1} (L_{N}(x) P_{N-1}(x) + Q_{N-1}(x)) dx.
= \int^1_{-1} Q_{N-1}(x) dx
\end{eqnarray}

如果将$Q_{N-1}(x)$也用勒让德多项式表达，有
\begin{equation}
Q_{N-1}(x) = \sum^{N-1}_{i=0} b_i L_i(x),
\end{equation}
而 $L_0(x)=1$，所以有
\begin{eqnarray}
\int^1_{-1} f(x)dx \approx \int^1_{-1} Q_{N-1}(x) dx　＝　\int^1_{-1} L_0(x) Q_{N-1}(x) dx
= 2 b_0,
\end{eqnarray}
所以只要我们知道了$b_0$的值，就得到了积分的近似值。
\end{frame}

\begin{frame}{高斯积分}
怎样求取$b_0$？我们能计算的是$f(x)$。。。
\begin{equation}
f(x) \approx L_{N}(x) P_{N-1}(x)　+ Q_{N-1}(x),
\end{equation}
$L_N(x)$在$(-1,1)$上有$N$个零点,我们把它们叫做$x_0, \cdots, x_{N-1}$，
则有
\begin{eqnarray}
f(x_i) =　L_{N}(x_i) P_{N-1}(x_i)　+ Q_{N-1}(x_i) = Q_{N-1}(x_i) = \sum^{N-1}_{j=0} b_j L_j(x_i),
\end{eqnarray}
$L_j(x_i), i,j=0,\cdots,N-1$构成一个方阵$L$，可以证明其行列式不为0，所以
\begin{equation}
b_0 = \sum^{N-1}_{i=0} f(x_i) L^{-1}_{i0},
\end{equation}
积分近似为
\begin{equation}
\int^1_{-1} f(x)dx \approx 2 \sum^{N-1}_{i=0} f(x_i) L^{-1}_{i0}
= \sum^{N-1}_{i=0} \omega_i f(x_i).
\end{equation}
\end{frame}

\begin{frame}{高斯积分}
\begin{equation}
\int^1_{-1} f(x)dx \approx 2 \sum^{N-1}_{i=0} f(x_i) L^{-1}_{i0}
= \sum^{N-1}_{i=0} \omega_i f(x_i).
\end{equation}
$x_0,\cdots,x_{N-1}$是固定的，$2L^{-1}_{i0}$也是固定的，所以我们可以求得它们，对任何$f(x)$应用。
下表是$N-1=5$的例子
\begin{table} [h]
\centering
\begin{tabular}{ccc}
\hline \hline
i & $x_i$ & $\omega_i = 2 L^-1_{i0}$ \\ 
\hline 
0 & -0.90617 98459 38664 & 0.23692 68850 56189  \\  
1 & -0.53846 93101 05683  & 0.47862 86704 99366  \\  
2 & 0.00000 00000 00000 & 0.56888 88888 88889  \\  
3 & 0.53846 93101 05683  & 0.47862 86704 99366  \\  
4 & 0.90617 98459 38664 & 0.23692 68850 56189  \\ 
\hline \hline
\end{tabular} 
\label{table:GaussLegendre5}
\caption{高斯勒让德积分：$L_5(x)$的零点，与相应的高斯积分权重。}
\end{table}
\end{frame}

\begin{frame}{高斯积分}
当然，积分区域不一定是$[-1,1]$，而是$[a,b]$，我们可以做如下变换
\begin{equation}
t = \frac{ x- (a+b)/2 }{ (b-a)/2 },
\end{equation}
则t的积分区域为$[-1,1]$,所以
\begin{eqnarray}
\int^b_{a} f(x)dx &=& \frac{b-a}{2} \int^1_{-1} f( \frac{a+b}{2} + \frac{b-a}{2}t ) dt 
\nonumber\\
&\approx& \frac{b-a}{2} \sum^{N-1}_{i=0} \omega_i f( \frac{a+b}{2} + \frac{b-a}{2}t_i ),
\end{eqnarray}
其中$t_i$为$L_N(x)$的零点，$\omega_i$为相应的权重，$N=5$时，如表所示。%\ref{table:GaussLegendre5}所示。
\end{frame}

\begin{frame}[fragile]{高斯积分}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<iomanip>

const double x_node[5]={
        -0.906179845938664,
        -0.538469310105683,
        0,
        0.538469310105683,
        0.906179845938664};
const double omega[5]={
        0.236926885056189,
        0.478628670499366,
        0.568888888888889,
        0.478628670499366,
        0.236926885056189};
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]
\lstset{firstnumber = last}
\begin{lstlisting}
double GaussLegendre5(double a, double b, double (*func)(double)){
        double x, y=0;
        for(int i=0;i<5;i++){
                x = x_node[i] * (b-a)/2 + (a+b)/2 ;
                y += omega[i] * func(x);
        }
        y *= (b-a)/2;
        return y;
}

double f(double x){
        return x*x*x*x*x;
}

int main(){

        double a=0;
        double b=1;
        cout<<"Integration of f(x) from "<<a<<" to "<<b<<" is ";
        cout<<setprecision(15)<<GaussLegendre5(a,b,f)<<endl;
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{高斯积分}
根据我们的理论分析，用高斯积分方法，相当于用$2N-1$阶多项式拟合被积函数，
\begin{eqnarray}
P_{2N-1}(x) = L_{N}(x) P_{N-1}(x) + Q_{N-1}(x),
\end{eqnarray}
在积分公式中，涉及到$L_N (x)$的$N$个零点，以及与之相关的$N$个权重。
\begin{eqnarray}
\int^b_{a} f(x)dx 
&\approx& \frac{b-a}{2} \sum^{N-1}_{i=0} \omega_i f( \frac{a+b}{2} + \frac{b-a}{2}t_i ),
\end{eqnarray}
上面我们写了$N=5$的代码，那么也就是说，我们用$9$阶多项式拟合被积函数。

那么我们得到一个结论：如果被积函数是多项式，阶数不大于9，我们的代码会非常精确；如果阶数大于9，理论上就有误差。

\end{frame}

\begin{frame}[fragile]{高斯积分}

\begin{table}
\begin{tabular}{ccc}
\hline \hline
$f(x)$ & $\int^1_0 f(x)dx$数值解 & 误差 \\ 
\hline 
$x$     & 0.5  & -2.22044604925031e-16 \\ 
 
$x^2$ & 0.333333333333333 & -1.66533453693773e-16 \\ 
 
$x^3$ & 0.25 & -1.66533453693773e-16 \\ 
 
$x^4$ & 0.2 & -1.66533453693773e-16 \\ 
 
$x^5$ & 0.166666666666667 & -1.38777878078145e-16 \\ 
 
$x^6$ & 0.142857142857143 & -1.11022302462516e-16 \\ 
 
$x^7$ & 0.125 & -1.2490009027033e-16 \\ 
 
$x^8$ & 0.111111111111111 & -1.11022302462516e-16 \\ 
 
$x^9$ & 0.0999999999999999 & -1.11022302462516e-16 \\ 
 
$x^{10}$ & 0.0909076593600402 & -1.43154905070386e-06 \\ 

$x^{11}$ & 0.0833254598135549 & -7.87351977839246e-06 \\

$x^{12}$ & 0.0768984805856498 & -2.45963374271296e-05 \\

\hline\hline
\end{tabular} 
\caption{5阶高斯积分：$\int^1_0 f(x) dx$，被积函数为多项式时的数值误差。}
\end{table}
\end{frame}

\begin{frame}[fragile]{高斯积分的零点与权重}
下面的函数改写自《Numerical Recipes in C》上的一段代码，根据积分区间$[a,b]$与阶数$n$，计算$[a,b]$区间内的节点，存放在$x[]$中，相应的权重存放在$w[]$中。
整段积分的值为
\begin{equation}
I_f \approx \sum^n_{i=1} w_i f(x_i).
\end{equation}

\begin{lstlisting}
/*
 * Given [a,b], returns abscissas x[] and weights w[], based on zeros of P_n(x)
 * Modified from Numerical Recipes
 */
void GauLeg(double a, double b, int n, double *x, double *w){

        double eps=1E-14;
        int m,j,i;
        double z1,z,xm,xl,pp,p3,p2,p1;

\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{高斯积分的零点与权重}
\begin{lstlisting}
        m=(n+1)/2;
        xm=0.5*(a+b);
        xl=0.5*(b-a);
        for(i=0;i<m;i++){
                z=cos(M_PI*(i+0.75)/(n+0.5));//starting guess
                do{                     //Newton's method to get the ith zero
                        p1=1.0;         // P_0
                        p2=0.0;         
                        for(j=0;j<n;j++){
                                p3=p2;  // P_{j-1}
                                p2=p1;  // P_{j}
                                p1= z*p2 + (z*p2-p3)*j/(j+1); 
                                // P_{j+1},based on recurrence relation:
                                // (j+1)P_{j+1} = (2j+1)P_j - j P_{j-1} 
                        }
                        pp=n*(z*p1-p2)/(z*z-1); // P'n(x_i), based on
                                    // (1-x^2)P'_n(x) = n P_{n-1}(x) - nxP_n (x)
                        z1=z;
                        z=z1-p1/pp; // Newton's method: s = x - f(x)/f'(x)
                }while(fabs(z-z1)>eps);
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{高斯积分的零点与权重}
\begin{lstlisting}
                x[i]=xm-xl*z;
                x[n-1-i]=xm+xl*z;
                w[i] = 2.0*xl/((1.0-z*z)*pp*pp); // formula of weights:
                                    // w[i] = 2/(1-x_j^2)/P'n(x_j)^2 for [-1,1]
                                    // multiplied by xl, for [a,b]
                w[n-1-i]=w[i];
        }
}
\end{lstlisting}
\end{frame}

\begin{frame}{高斯积分}
在实际应用中——我们再重复一遍——对于震荡比较厉害的被积函数，Newton-Cortes积分更加安全。\\
\vspace{0.5cm}
如果只取有限的几个点，进行高斯积分，是危险的。
这是因为，高斯积分的出发点是用高阶多项式拟合被积函数，而震荡很厉害的被积函数，可能需要非常高阶的多项式才能拟合好。。。\\
\vspace{1cm}
对于比较平缓的被积函数，用多项式拟合往往效果不错，所以可以使用高斯积分，只需取很少的几个点，可能非常高效。\\
\vspace{1cm}
当然，你也可以将积分区间分割成很多小份，然后在每个小区间里进行高斯积分。
\end{frame}

\begin{frame}
\frametitle{目录}
	\begin{itemize}
		\item { 第四章：数值求导与多项式拟合 }
		\vspace{1cm}
		\item 第五章：数值积分
		\vspace{1cm}
		\item {\color{blue}第六章：方程求根 }
        \vspace{1cm}
        \item 第七章：寻找函数极值
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{方程求根}
方程求根，即求解 
\begin{equation}
f(x)=0,
\end{equation}
实际计算中经常遇到没有解析解的方程，只能求数值解，比如
\begin{eqnarray}
\cos x -x =0,
\end{eqnarray}
只能求数值解。
有时$f(x)$甚至没有解析形式，可能是一系列计算的结果。
\end{frame}

\begin{frame}[fragile]{找所有的根？找一个根？}

\begin{itemize}
\item 方程有多少个根？

\begin{itemize}
\item 理论上判断
\item 等步长取点判断趋势(画个图看看)
\end{itemize}

\item 找某一个根？

\begin{itemize}
\item 二分法
\item 牛顿法(切线法)
\item 割线法
\end{itemize}

\end{itemize}

\end{frame}

\begin{frame}[fragile]{二分法}	
已知区间 $(a,b)$ 内有 $f(x)=0$ 的根，且 $f(a)*f(b)<0$，可以如下逼近方程的根。
\begin{itemize}
\item [1] 取$c=(a+b)/2$，计算$f(c)$。
\item [2] 若$f(c)f(a) \leq 0$，说明$[a,c]$ 区间内存在根。将$c$的值赋给$b$，回到步骤１。
\item [3] 否则，即$f(c)f(a) > 0$，说明$f(c)*f(b)<0$，区间$[c,b]$内存在根，将$c$的值赋给$a$，回到步骤1。
\end{itemize}
每次迭代，区间都缩小一半。如此循环，不断缩小“嫌疑区间”，就可以逼近根的真实值。
当区间小于目标精度时，即可取区间的中点，作为方程的根。
\end{frame}
%由于每次循环中区间的长度是上一次循环的一半，所以收敛速度是线性的。

\begin{frame}[fragile]{二分法}	
\begin{lstlisting}
#include<iostream>
#include<cmath>
#include<fstream>
#include<string>
using namespace std;
#define max_iteration 1E6
/*
 * findroot: 在[a,b]区间内寻找 func(x)=0　的根，并返回根的值
 * 算法：二分法
 */
double bisection(double (*func)(double), double a, double b, double precision){

        double ya=func(a), yb=func(b);
        double c=(a+b)/2, yc;
        int i;

        if( ya* yb > 0){//如果 func(a) * func(b) >0, 则[a,b]区间内不一定有解，拒绝任务
                cout<<"error: invalid input [a,b] for findroot: f(a)*f(b)>0 "<<endl;
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{二分法}	
\lstset{firstnumber=last}
\begin{lstlisting}
        for(i=0;i<max_iteration;i++){

                c = (a+b)/2;//区间中点
                yc = func(c);

                if( yc * ya <= 0 ){//如果 yc*ya<=0，说明[a,c]中有一个根
                        b=c;
                }
                else{//如果 yc*ya>0，说明[c,b]中有一个根
                        a=c;
                }
                if( fabs(a-b) < precision ){
                        return (a+b)/2;
                }
        }
        cout<<"bisection method: failed to find a root, after "<<max_iteration<<" iterations.\n";
        exit(1);
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{二分法}	
\lstset{firstnumber=last}
\begin{lstlisting}
double g(double x){

        return x*x - 4*x -5;
}

int main(){

        double a=0, b=1, c;

        while( g(a) * g(b) >0 ){
                cout<<"Input valid boundaries inside which we get a root."<<endl;
                cout<<"a=";
                cin>>a;
                cout<<"b=";
                cin>>b;
        }

        cout<<" The root is x="<<(c=bisection(g,a,b,1E-9))<<" f(x)="<<g(c)<<endl;
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{牛顿法}	
如果方程的根是$s$，$f(s)$可以对任意$x$处做泰勒展开
\begin{equation}
0 = f(s) = f(x) + f'(x)(s-x) + \frac{ f''(x) }{ 2 } (s-x)^2 + \cdots
\end{equation}
当$s$很接近$x$时，忽略高阶项，有
\begin{equation}
s \approx x - f(x)/f'(x),
\end{equation}
受这个启发，可以构造一个迭代式
\begin{eqnarray}
x_{n+1} = x_n - \frac{ f(x_n) }{ f' (x_n)}
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{牛顿法}	
这个迭代具有几何的意义。
在$y=f(x)$曲线上$(x_n, f(x_n))$处作切线，与$x$轴的交点正是上面定义的$x_{n+1}$。
所以牛顿饭又叫做切线法。
迭代下去，当$f(x_{n+1})$小于某个很小的数时，即可把$x_{n+1}$当做方程根的近似值。
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{NewtonMethodDemon.eps}
\caption{牛顿法逼近方程的根：切线与x轴交点为下一个迭代值。}
\label{figure:NewtonMethodDemon}
\end{figure}
\end{frame}

\begin{frame}[fragile]{收敛速度}	
若真实值为 $\alpha$，$x_n$ 的误差为 $\epsilon_n = x_n - \alpha$，则
\begin{eqnarray}
x_{n+1} &=& x_n - \frac{ f(x_n) }{ f' (x_n)} \nonumber\\
				  &\approx& \alpha + \epsilon_n - \frac{ \epsilon_n f'(\alpha) + \frac{1}{2} \epsilon^2_n f''(\alpha)}{f'(\alpha) + \epsilon_n f''(\alpha)} \nonumber\\
				  	&\approx& \alpha + \frac{1}{2} \frac{ f''(\alpha) }{ f'(\alpha)} \epsilon^2_n,
\end{eqnarray}
即
\begin{eqnarray}
\epsilon_{n+1} \approx \frac{1}{2} \frac{ f''(\alpha) }{ f'(\alpha)} \epsilon^2_n,
\end{eqnarray}
所以，如果$\epsilon_n$是$10^{-3}$，$\frac{1}{2} \frac{ f''(\alpha) }{ f'(\alpha)} \approx 1$的话，$\epsilon_{n+1}$约为$10^{-6}$量级。
一般来说，牛顿法收敛得比二分法快。
\end{frame}

\begin{frame}[fragile]{收敛速度}
但是牛顿法不一定收敛，一个典型的例子是，$f'(x_n) = 0$，迭代公式就没法用了。
可以先用二分法找到一个比较小的区间，再用牛顿法求方程的根。

另外，牛顿法需要使用一阶导数值，如果$f'(x)$没有解析值，而需要通过差分，可能引入新的误差，这种情况下，最好不要用牛顿法，可以使用后面介绍的割线法。
\end{frame}

\begin{frame}[fragile]{牛顿法}
\begin{lstlisting}
/*
 * Newton's method: given a function (*func), which gives value and also the 1st derivative,
 * a suspicious zone [a,b], it finds a root in the zone within precision, and returns the root.
 */
double Newton_findroot(void (*func)(double, double &, double &), double a, double b, double precision){

        double x, dx, y, dy;
        int i;

        x = (a+b)/2;
        
        for(i=0;i<max_iteration;i++){
                func(x, y, dy);
                if( fabs(dy)<1E-9 ){
                        cout<<" error: f'=0 in Newton's root-finding method. \n";
                        exit(1);
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{牛顿法}
\lstset{firstnumber=last}
\begin{lstlisting}
                dx = y/dy;
                x -= dx;
                if( (x-a)*(x-b) >0 ){
                        cout<<" error: x gets out of the given zone [a,b] in Newton's method.\n";
                        exit(1);
                }
                if( fabs(dx) < precision ){
                        return x;
                }
        }
        cout<<" Newton's method: after "<<max_iteration<<" steps, failed to find a root.\n";
        exit(1);
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{形式参数、变量别名}
在前面的代码中，我们使用了变量别名，所以做如下说明。
\begin{lstlisting}
void g(double x){
    x=x+1;
}
int main(){
    double a=1;
    g(a);                //调用g函数
    cout<<"a="<<a<<endl; //输出a的值
    return 0;
}
\end{lstlisting}
虽然g函数中x的值改变了，但是主函数中仍会输出“a=1”。

g的定义中x是形式参数，函数被调用时，它开辟一个新的内存单元，叫做x，然后进行操作，函数调用结束以后释放掉这块内存。

在主函数中调用g(a)时，a是实际参数，a的值被传给了新的内存单元x，g函数调用结束以后，a的值并没有被改变，只是值被读取了。
\end{frame}

\begin{frame}[fragile]{形式参数、变量别名}
\begin{lstlisting}
void g(double & x){//变量别名
    x=x+1;
}
int main(){
    double a=1;
    g(a);                //调用g函数
    cout<<"a="<<a<<endl; //输出a的值
    return 0;
}
\end{lstlisting}
在g的参数中，如果使用double \& x，则表示x是传入的a的别名。在g中对x的操作，就是对a的操作，所以调用g(a)以后a的值被改变了。
\end{frame}

\begin{frame}[fragile]{割线法}
如果$f'(x)$很难得到，则可以用割线代替牛顿法中的切线。
\begin{itemize}
\item [1] 取两个初始点$x_0, x_1$，通过$(x_0, f(x_0)), (x_1, f(x_1))$两点作直线，交于 $x$ 轴，记该点为 $x_2$。
\item [2] 取点$x_1, x_2$，通过$(x_1, f(x_1)), (x_2, f(x_2))$两点作直线，交于 $x$ 轴，记该点为 $x_3$。
\item [3] $\cdots$
\end{itemize}
迭代公式为
\begin{eqnarray}
x_{n+1} = x_n - \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} f(x_n).
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{割线法}
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{SecantMethodDemon.eps}
\caption{割线法逼近方程的根：割线与x轴交点为下一个迭代值。}
\label{figure:SecantMethodDemon}
\end{figure}
\end{frame}

\begin{frame}[fragile]{割线法}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

#define max_iteration 1E6

double g(double x){
        return x*x-1;
}

/*
 * SecantMethod: given function func, suspicious zone [a,b], returns a root of func(x)==0
 */
double SecantMethod(double (*func)(double), double a, double b, double precision ){

        double c, fc, d, fd, e, fe;
        int i;

        c = (a+b)/2;//set starting points
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]
\lstset{firstnumber=last}
\begin{lstlisting}
        fc = func(c);
        d = (c+b)/2;
        fd = func(d);
        for(i=0;i<max_iteration;i++){
                e = d - (d-c)/(fd-fc)*fd;
                if( (e-a)*(e-b) >0 ){
                        cout<<"SecantMethod: error! iteration gets out of the given zone [a,b].\n";
                        exit(1);
                }
                fe = func(e);
                if( fabs(e-d)<precision && fabs(fe)<precision ){
                        return e;
                }
                c=d;
                fc=fd;
                d=e;
                fd=fe;
        }
        cout<<"After "<<max_iteration<<" iterations, failed to find one root"<<endl;
        exit(1);
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{割线法}
\lstset{firstnumber=last}
\begin{lstlisting}
int main(){

        double a=0, b=10;
        cout<<" root: x="<<SecantMethod(g, a, b, 1E-6)<<endl;
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{收敛速度}
记方程的根的真实值为 $\alpha$，$x_{n-1}, x_n$ 的误差为 $\epsilon_{n-1} = x_{n-1} - \alpha, \epsilon_n = x_n - \alpha$，则
\begin{eqnarray}
x_{n+1} &=& \alpha + \epsilon_n - \frac{ \epsilon_n - \epsilon_{n-1} } { f'(\alpha) (\epsilon_n - \epsilon_{n-1}) + \frac{ f''(\alpha) }{2} (\epsilon^2_n - \epsilon^2_{n-1})}
\nonumber\\
&\approx& \alpha + \frac{ f''(\alpha) }{2f'(\alpha)} \epsilon_n \epsilon_{n-1},
\end{eqnarray}
即
\begin{equation}
\epsilon_{n+1} \approx \frac{ f''(\alpha) }{ 2f'(\alpha) } \epsilon_n \epsilon_{n-1}.
\end{equation}
\end{frame}

\begin{frame}[fragile]{收敛速度}
\begin{equation}
\epsilon_{n+1} \approx \frac{ f''(\alpha) }{ 2f'(\alpha) } \epsilon_n \epsilon_{n-1}.
\end{equation}
假设$ \epsilon_1 = \epsilon_0^r $，则
\begin{eqnarray}
\epsilon_0 &=& \epsilon_0^{1} \nonumber\\
\epsilon_1 &\propto& \epsilon_0^{r} \nonumber\\
\epsilon_2 &\propto& \epsilon_0^{r +1} \nonumber\\
\epsilon_3 &\propto& \epsilon_0^{2r +1} \nonumber\\
\epsilon_4 &\propto& \epsilon_0^{3r +2} \nonumber\\
\epsilon_5 &\propto& \epsilon_0^{5r +3} \nonumber\\
\epsilon_6 &\propto& \epsilon_0^{8r +5} \nonumber\\
\epsilon_7 &\propto& \epsilon_0^{13r +1} \nonumber\\
\cdots
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{收敛速度}
其中出现了两列裴波那契数列 $1, 1, 2, 3, 5, 8, \cdots$，
显然有
\begin{eqnarray}
\epsilon_n \propto \epsilon_0^{a_{n-1} r + a_{n-2}}.
\end{eqnarray}
而裴波那契数列的通项公式为 
\begin{eqnarray}
a_n = \frac{1}{\sqrt{5}} \left[ (\frac{1+\sqrt{5}}{2})^n - (\frac{1-\sqrt{5}}{2})^n \right],
\end{eqnarray}
且
\begin{eqnarray}
\lim\limits_{n \rightarrow \infty} \frac{ a_n} {a_{n-1}} = \frac{ 1+ \sqrt{5}}{2} \approx 1.618033,
\end{eqnarray}
所以 $n \rightarrow \infty$ 时，
\begin{eqnarray}
\epsilon_{n+1} \propto \epsilon_{n}^{1.618033}.
\end{eqnarray}
显然，割线法没有牛顿法收敛得快。
\end{frame}

\begin{frame}[fragile]{Stewenson方法}
\begin{eqnarray}
x_{n+1} = x_n - \frac{ f(x_n)^2 }{f(x_n) - f(x_n - f(x_n))}.
\end{eqnarray}
即用 $(x_n, f(x_n)), (x_n - f(x_n), f(x_n - f(x_n)))$ 两点作割线，与割线法类似。
但是 Stewenson 方法收敛更快，和牛顿法一样，都是2阶收敛。
\end{frame}

\begin{frame}[fragile]{Stewenson方法}
假设方程的根的精确值为 $\alpha$，$x_n$ 的误差为 $ \epsilon_n = x_n - \alpha$，
经过推导，可以证明
%\begin{eqnarray}
%&& x_{n+1} \approx \alpha + \epsilon_n \nonumber\\
%&& - \frac{ (f'(\alpha) \epsilon_n + \frac{1}{2} f''(\alpha) \epsilon^2_n)^2 }
%{ f'(\alpha) \epsilon_n + \frac{1}{2}f''(\alpha)\epsilon^2_n - \left[ f'(\alpha)(\epsilon_n -f'(\alpha)\epsilon_n -\frac{1}{2}f''(\alpha)\epsilon^2_n) 																																	+ \frac{1}{2}f''(\alpha)(\epsilon_n - f'(\alpha)\epsilon_n -\frac{1}{2}f''(\alpha)\epsilon^2_n)^2 \right]}.																																	\nonumber\\
%\end{eqnarray}
%上式最后一项分子取到 $\epsilon^3_n$ 项，分母取到 $\epsilon^2_n$ 项，得到
%\begin{eqnarray}
%x_{n+1}  \approx \alpha + \frac{ f''(\alpha) (1- f'(\alpha)^2)}{f'(\alpha)^2} \epsilon^2_n,
%\end{eqnarray}
%即
%
\begin{eqnarray}
\epsilon_{n+1} \approx  \frac{ f''(\alpha) (1- f'(\alpha)^2)}{f'(\alpha)^2} \epsilon^2_n.
\end{eqnarray}
即Stewenson方法是二阶收敛的。
所以它不需要用到导数，但也可以和牛顿法收敛得差不多快。
缺点是公式稍微复杂一点点。
\end{frame}

\begin{frame}{应用：氘核}
氘原子核由1个质子和1个中子组成，我们可以将它们之间的相互作用简化为有限深势井，即
\begin{equation}
V(r_{pn}) = \left\{
\begin{aligned}
&-V_0, ~~ &r < a, \\
&0, ~~ &else,
\end{aligned}
\right.
\end{equation}
其中$r$是质子中子之间的距离，$a$为$fm$数量级。
两体系统可以由质心坐标、相对坐标来描述。
质心坐标描述了原子核的平动，比如氘原子在水分子中，随着水分子在空气中做热运动。
相对坐标则描述了原子核内部的结构：e.g.基态时质子中子的相对距离。

相对坐标的薛定谔方程可以写作
\begin{equation}
(-\frac{\hbar^2}{2\mu} \nabla^2 + V(r) ) \Psi = E \Psi,
\end{equation}
其中$\mu　\equiv \frac{ m_p m_n }{m_p + m_n}$为折合质量。
\end{frame}

\begin{frame}{氘核}
分离变量以后，可以得到径向方程
\begin{equation}
(-\frac{\hbar^2}{2\mu} \frac{ d^2 }{d r^2} + V(r) + \frac{\hbar^2}{2\mu} \frac{l(l+1)}{r^2} ) u(r) = E u(r),
\end{equation}
而波函数为$\Psi = \frac{u(r)}{r} Y_{lm}(\theta, \phi)$.

在我们关于$V(r)$的假定下，体系的基态为$l=0$的束缚态。$l=0$时，径向方程为
\begin{equation}
(-\frac{\hbar^2}{2\mu} \frac{ d^2 }{d r^2} + V(r)) u(r) = E u(r),
\end{equation}
在$r <a$时，上述方程为
\begin{equation}
\frac{ d^2 }{d r^2}u(r)　+ K^2 u(r) = 0,
\end{equation}
其中$K = \frac{ \sqrt{ 2 \mu (E+V_0) } }{\hbar}$。
\end{frame}

\begin{frame}{氘核}
在$r <a$时，
\begin{equation}
\frac{ d^2 }{d r^2}u(r)　+ K^2 u(r) = 0,
\end{equation}
其中$K = \frac{ \sqrt{ 2 \mu (E+V_0) } }{\hbar}$。
$r > a$时，$u(r)$满足方程
\begin{equation}
\frac{ d^2 }{d r^2}u(r)　- \kappa^2 u(r) = 0,
\end{equation}
其中$\kappa = \frac{ \sqrt{ - 2 \mu E } }{ \hbar}$。

所以
\begin{equation}
u(r) = \left\{　
\begin{aligned}
& A\sin (Kr) + B \cos (Kr), & 0<r\leq a, \\
& Ce^{-\kappa r }, & r>a
\end{aligned}
\right.
\end{equation}
\end{frame}

\begin{frame}{氘核}
由于$\Psi = \frac{u(r)}{r} Y_{lm}(\theta, \phi)$在$r\rightarrow0$时为有限值，在$r=a$处连续，且一阶导数连续，得到
\begin{eqnarray}
B &=& 0, \\
A\sin(Ka) &=& C e^{-\kappa a}, \\ 
AK\cos(Ka) &=& -C \kappa e^{-\kappa a}.
\end{eqnarray}
得到
\begin{equation}
K \cot(Ka) = - \kappa,
\end{equation}
即本证能量为下面的函数的零点，
\begin{equation}
f(E) = \sqrt{2\mu(E+V_0)} cot( \sqrt{2\mu a^2 (E+V_0) } / \hbar ) + \sqrt{-2\mu E}. 
\end{equation}
\end{frame}

\begin{frame}{氘核}
取自然单位制$\hbar = c = 1$, $fm \approx 5 GeV^{-1}$，折合质量为$\mu = 0.5GeV$，
\begin{equation}
f(E)　(MeV) = \sqrt{1000(E+V_0)} cot( \sqrt{0.025 a^2 (E+V_0) } ) + \sqrt{-1000 E}, 
\label{eqn:deuteron-f(E)}
\end{equation}
其中$a$的单位为$fm$，$E, V_0$的单位为$MeV$。

如果取$a=2fm$，$V_0=35MeV$，$f(E)$的图如下
\begin{figure}
\centering
\includegraphics[width=0.5\linewidth]{deuteron}
%\caption{Deuteron s-wave bound state: $f(E)$ in (\ref{eqn:deuteron-f(E)}) v.s. E in MeV.}
\label{fig:deuteron}
\end{figure}
对应的是氘核的结合能，约2.2MeV.

\end{frame}

\begin{frame}
\frametitle{目录}
	\begin{itemize}
		\item { 第四章：数值求导与多项式拟合 }
		\vspace{1cm}
		\item 第五章：数值积分
		\vspace{1cm}
		\item {第六章：方程求根 }
        \vspace{1cm}
        \item {\color{blue}第七章：寻找函数极值 }
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{寻找函数极值}
\begin{itemize}
\item 有时候一个物理量是一个参数的函数$f(x)$，而我们需要找到最优参数，使得$f$取得最小值（或者最大值）。
\item 比如，量子力学中的变分法，就是一个求极小值的问题。
\item 最简单的方法是，我们在参数的取值范围内取一系列离散的点$x_i$，取$f(x_i)$中的最小值。
\end{itemize}

\begin{figure}[h]
\centering
\includegraphics[width=0.6\linewidth]{optimization_demon}
\caption{示意图：寻找函数极值，理论上的极值为x=-0.5。}
\label{fig:optimization_demon}
\end{figure}
\end{frame}

\begin{frame}[fragile]{寻找函数极值}
如果希望做得更加精确一点，可以用算法来逼近极小值。
\begin{shaded*}
\begin{itemize}
\item 第一步是找到一个区间，确保这个区间内有极小值。
\item 第二步是不断缩小这个区间，逼近极小值。
\end{itemize}
\end{shaded*}

\begin{figure}[h]
\centering
\includegraphics[width=0.6\linewidth]{optimization2}
%\caption{示意图：寻找函数极值，理论上的极值为x=-0.5。}
\label{fig:optimization2}
\end{figure}

\end{frame}

\begin{frame}[fragile]{寻找函数极值}
第一步比较容易，可以如下实现
\begin{shaded*}
\begin{itemize}
\item [1] 给定一个初始区间$[a,b]$, 比较$f(a), f(b)$两个值，如果$f(a)<f(b)$，调换$a,b$的值。寻找极小值的方向为$ a \rightarrow b \rightarrow $。
\item [2] 定义$c = b + w*(b-a)$, 其中$w$为一个常系数，代表了每次前进的相对步长大小，常取作 1.618。
\item [3] 如果$ f(c) < f(b) $，说明函数值还在减小，则令$a=b, b=c$，回到步骤2。
\item [4] 如果$ f(c) > f(b) $，说明$ [a,c] $区域内有极小值。
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{寻找函数极值}
\begin{figure}[h]
\centering
\includegraphics[width=0.7\linewidth]{optimization3}
\caption{示意图：寻找函数极值，理论上的极值为x=-0.5。}
\label{fig:optimization3}
\end{figure}
比如，在上图中，若初始区间为$[a,b]=[0,0.1]$，则会在第1步得到$[a,b]=[0.1,0]$，在第2步得到$c=-0.1618$，反复执行几次2,3步骤以后，在第4步得到一个区间，确认该区间内有极小值。
\end{frame}

\begin{frame}[fragile]{寻找函数极值}
第二步是缩小“嫌疑空间”，逼近极小值。
如果初始区间为 $[a,c]$，其中有一点 $b$，满足 $f(a)>f(b), f(c)>f(b) $，可以如下实现完成第二步
\begin{shaded*}
\begin{itemize}
\item [1] 取 $[a,c]$ 区间内另一点 $d$，一般取 $[a,c]$ 的某个比例分割点，比如黄金分割点$d=a+0.618(c-a)$，求得 $f(d)$ 的值。
\item [2] 如果$d \in (b,c)$：①如果 $f(d)<f(b)$，则 $a=b, b=d$；② $f(d)>f(b)$，则 $c=d$。
\item [3] 如果$d \in (a,b)$：①如果 $f(d)<f(b)$，则 $c=b, b=d$；② $f(d)>f(b)$，则 $a=d$。
\item [4] 如果$|a-c|$小于某个指定的精度，就将$(a+c)/2$当作极小值位置。
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{寻找函数极值}
实现上述功能的代码如下
\begin{lstlisting}
#include<iostream>
using namespace std;

#include<cmath>

#define max_iteration 1E4

/*
 * findmin: 给定函数func, 初始区间[a,b]，在精度要求precision下，寻找一个极小值点，并返回该点坐标
 */
double findmin(double (*func)(double), double a, double b, double precision){

        double temp,c,d;
        int i;
        if( func(a)<func(b) ){//寻找极小值的方向为：a->b->
                temp=a;
                a=b;
                b=temp;
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{寻找函数极值}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(i=0;i<max_iteration;i++){//寻找嫌疑区间
                c= b + 1.618*(b-a);
                if( func(c) > func(b) )break;
                a=b;
                b=c;
        }
        if(i==max_iteration){//10000次迭代以后，没有找到嫌疑区间，报错
                cout<<"after "<<max_iteration<<" iterations, didn't find a zone for fmin"<<endl;
                exit(1);
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{寻找函数极值}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(i=0;i<max_iteration;i++){
                if( fabs(a-c)<precision )break;
                d=0.38197*a + 0.61803*c;
                if( (d-a)*(d-b)<0){//如果d在a,b之间
                        if( func(d) < func(b) ){
                                c=b;
                                b=d;
                        }
                        else{
                                a=d;
                        }
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{寻找函数极值}
\lstset{firstnumber = last}
\begin{lstlisting}
                else if( (d-a)*(d-b)>0 ){//如果d在b,c区间
                        if( func(d) < func(b) ){
                                a=b;
                                b=d;
                        }
                        else{
                                c=d;
                        }
                }
                //下面四行帮忙查错：观察每次迭代的结果
                //cout<<"-----------------------------------"<<endl;
                //cout<<"a="<<a<<", f(a)="<<func(a)<<endl;
                //cout<<"b="<<a<<", f(b)="<<func(b)<<endl;
                //cout<<"c="<<a<<", f(c)="<<func(c)<<endl;
        }
        return (a+c)/2;
}
\end{lstlisting}
\end{frame}

\begin{frame}{作业}
１．利用Richardson算法，计算$e^x$在$x=10$处的二阶导数，计算$D^0_3$，并计算误差。

*2. 牛顿也写了一个多项式插值公式，他定义
\begin{eqnarray}
P_n (x) &=& a_0 + a_1 (x-x_0) + a_2(x-x_0)(x-x_1) \nonumber\\
        && + \cdots + a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1}),
\end{eqnarray}
它的前1项通过第一个数据点，前2项通过前两个数据点，...，前n+1项在一起
通过前n+1个数据点 $(x_0,y_0), (x_1,y_1), \cdots, (x_n,y_n)$。

1)试证明
\begin{equation}
a_n = \frac{ y_n - P_{n-1} (x_n) }{ (x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1}) },
\end{equation}

2)试用递归函数编写一个程序，自动读入数据，用牛顿插值公式进行拟合，并给出指定x值处的插值$P_n(x)$。
\end{frame}

\begin{frame}{作业}
3．粒子从$t=0s$到$t=10s$，速度为$v(t)=10t-t^2 (m/s)$，请用梯形法，积分得到粒子走过的路程，并与解析解相比较，误差为多大？

*4.　编写10阶高斯-勒让德积分，计算$\int^1_0 e^x dx$，并与理论值比较，得出误差。

5.　将$[0,\pi]$分割为10000份，用中值法、梯形法、辛普生方法分别积分$\int^\pi_0 \sin(x) dx$，比较误差，并进行简单分析。

*6.　试讨论：怎样构造自适应算法+梯形法进行积分？

7.试选取适当的区间，找到$f(x)=x(x-1)(x-2)(x-3)$的两个极小值。

*8.寻找$f(x)=x(x-1)(x-2)(x-3)$的一个极大值。

9.试编写辛普生积分的openmp并行程序，并自选一个任务，记录并行v.s.单线程完成任务所需时间。

10.任选一种方法，找到$f(x) = x^2 -4x -5$的两个零点，要求精度在$10^{-9}$以内，并给出迭代次数。

*11.试用Stewenson方法，找到$f(x) = x^2 -4x -5$的两个零点，要求精度在$10^{-9}$以内，并给出迭代次数。

\end{frame}




\end{document}